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Abstract 

The dynamics of a three-level atom in a cascade configuration with both transitions coupled to 
a single structured reservoir of quantized field modes is treated using Laplace transform methods 
applied to the coupled amplitude equations. Results are also obtained from master equations by 
two different approaches, that is, involving either pseudomodes or quasimodes. Two different types 
of reservoir are considered, namely a high-Q cavity and a photonic band-gap system, in which the 
respective reservoir structure functions involve Lorentzians. Non-resonant transitions are included 
in the model. In all cases non-Markovian behaviour for the atomic system can be found, such as 
oscillatory decay for the high-Q cavity case and population trapping for the photonic band-gap 
case. In the master equation approaches, the atomic system is augmented by a small number 
of pseudomodes or quasimodes, which in the quasimode approach themselves undergo Markovian 
relaxation into a flat reservoir of continuum quasimodes. Results from these methods are found to 
be identical to those from the Laplace transform method including two-photon excitation of the 
reservoir with both emitting sequences. This shows that complicated non-Markovian decays of an 
atomic system into structured EM field reservoirs can be described by Markovian models for the 
atomic system coupled to a small number of pseudomodes or quasimodes. 

PACS numbers: 42.50.Ct, 03.65.Yz 



1 



I. INTRODUCTION 



The study of open quantum systems, in which the interaction of the quantum system 
with the environment is taken into account ,Js of fundamental importance in several areas of 
physics. One such area is quantum optics 1], where the quantum system may be an atom 
(or an atom combined with a single mode laser field) and the environment with which the 
atom interacts may be a continuum of quantised modes of the electromagnetic fieldfor a set 
of vibrational modes in a solid). Another area is quantum measurement theory 0, ^, where 
the environment includes the macroscopic measuring apparatus, whose states "record" the 
results of measurements on the quantum system itself. Explaining the emergence of the 
classical world is related to the way the presence of the environment destroys coher- 

ence between certain states of macroscopic systems, resulting in only the "pointer states" 
remaining stable — any quantum superpositions of pointer states are rapidly converted into 
mixed states. A third area is that of degenerate quantum gases P], where the quantum 
system is an atomic Bose-Einstein condensate or atomic Fermi gas, and the environment 
may consist of those atoms in thermally excited states, as distinct from being in the macro- 
scopically occupied condensate state or in states within the Fermi surface. As in some of 
these examples, the quantum system itself may be microscopic or macroscopic. 

In many cases the interaction between the quantum system and the environment (also 
referred to as the bath or reservoir) involves coupling constants and reservoir mode densities 
that have a slowly varying frequency dependence. In such cases the reservoir structure 
function (which is the product of the mode density with the square of the coupling constant) 
is also slowly varying, resulting in the reservoir correlation time (which is the inverse of 
the bandwidth for the reservoir structure function) being very short compared to the time 
scales for dynamic evolution of the quantum system. This situation enables the dynamical 
behaviour of the quantum system and its interaction with the environment to be described 
by Markovian master equations (see e.g. Q). A large literature exists where the dynamical 
behaviour of quantum systems coupled to the environment has been successfully explained 
via Markovian master equations and related methods, and indeed much of quantum optics 
(see 111, 01) falls into this category. 

In recent years however, there has been an interest in open quantum systems where the 
conditions required for Markovian behaviour do not necessarily apply. Cases include atom 
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lasers quantum Brownian motion j8|, |9|, systems with conditioned evolution (such as 
associated with photodetection) [l^, decoherence in large scale quantum computers 11, 12, 
Q and in macroscopic systems generally Q, Q, Q. A further non-Markovian situation 
occurs for atomic systems coupled to structured reservoirs of electromagnetic (EM) field 
modes, where either the coupling constants or the mode density (or both) change rapidly 
with frequency. This situation can occur for atoms in high-Q cavities or in photonic band- 
gap (PBG) systems, and a general review of such situations is given by Lambropoulos et al. 



A number of methods for treating non-Markovian problems have been formulated. These 

3|; Heisenberg equations of motion l3: stochastic 
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tion operator master equation 

wave function methods for non-Markovian processes [8|, |10|, |2j, |25 , 
ods based on the essential states approximation or resolvent operators [1 
pseudomode app roach ^[:_ Fano diagonalisation Q, 1^^; and various short time scale 



methods 



41 




SC I; m eth- 

m 



the 



sense, but often simpler to appl y a nd interpret. 



In recent work 
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46|) we have studied non-Markovian processes for atomic 



systems coupled to a structured reservoirs of EM field modes, with applications to high- 
Q cavities and PBG systems Q|. Our approach has been based on the essential states 
approximation, Fano diagonalisation and pseudomodes. In Ref. the case of a three level 
cascade (or ladder) system coupled to a structured reservoir of EM modes was treated using 
the essential states approach, following a method similar to jssl. The general equations 
for the coupled amplitudes and their Laplace transforms were obtained, and the resulting 
integral equation solved via numerical methods, based on discretising frequency space to 
give a matrix equation that is equivalent to the original integral equation and also utilising 
analytical continuation in the complex s plane. The decay of the initially excited upper 
state was found as a function of time. Our application in [4^ was restricted to the case 
where the structured reservoir was due to a high-Q cavity. In the present paper, we apply 
the same method also to the case where the structured reservoir is due to a model photonic 
band-gap. In addition some further high-Q cavity cases are examined, where the cavity 
resonance is detuned from the atomic transition frequency. Other non-Markovian studies of 
the cascade system in a PBG system have been carried out in Refs. 0,13]. In these papers 
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a discontinuous mode density is used, whereas we represent the PBG by a difference of two 
Lorentzian functions. Quantitative comparisons between their results and ours are therefore 
not possible, only general features such as population trapping effects are in common. For 
completeness and to set out the notation, we include here in Section |n] the key integral 
equations for the cascade system decay, based on the essential states approach, that were 
derived in 4^. Results for both the high-Q cavity and PBG cases are presented in Section 
nil The dynamics is also treated via master equations by two different approaches, that 
is, involving either pseudomodes or quasimodes. In the master equation approach, the 
atomic system is augmented by a small number of pseudomodes or quasimodes, which in 
the quasimode approach themselves undergo Markovian relaxation into a fiat reservoir of 



continuum quasimodes.^ 
pseudomode approach 



he master equation theory is presented in Section IIV[ with the 
34L Isst being described in Subsection IIV Al and a treatment based 
on quasimodes \^ given in Subsection IIVBI Results from the master equation method 
are compared to those from the Laplace transform method. A summary of the paper is set 
out in Section El 



II. ESSENTIAL STATES THEORY 



A. Model System 

The system consists of a three level cascade atom, shown in figure ^ coupled to a bath 
(or reservoir) of bosonic modes. The bath modes would usually be associated with the 
quantum EM field, but other bosonic baths (such as lattice vibrations in a solid) might 
realise a similar model. For convenience, we will refer to the bath quanta as photons, but 
other quanta could be involved (such as phonons in the lattice vibration case). 

The atomic states are illustrated in figure Q the energy difference between states 1 and 
is hui, and between 2 and 1 the difference is ^ci;2. For later use we also define a frequency 
offset S, such that 

5=iLU,-UJ2)/2. (1) 

The offset determines the energy difference between the two transitions. If the offset is 
large, we would expect the system to behave as though there are two independent reservoirs. 
However, for relatively small offsets a photon may well be emitted from one transition, and 
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FIG. 1: Energy levels and energy differences in the three level system. The off-set S represents the 
difference between the energy of level 1 and the midpoint between levels and 2. 

later re-absorbed on the other transition. The bath modes have frequencies cux, raising and 
lowering operators a\ and a\ and the density of modes is px. Frequency dependent coupling 
constants g\i and gx2 specify the coupling between atomic transitions — 1 and 1 — 2 and 

the reservoir. 

The Hamiltonian for the system in the rotating wave approximation is given by: 



H = h 



a;i|l)(l| + (a;i +a;2)|2)(2| + ^a;A4aA 



+ E [^Ai (4|0)(1| + aA|l)(0|) + gx2 (al|l)(2| + aA|2)(l|) 



(2) 



Initially the atom is assumed to be in the upper state |2) and the bath modes are empty of 
photons, so the initial state vector is 



I^(0)) = |2)|0a). 



(3) 



The atom plus bath system evolves according to the Schrodinger equation. The initial state 
assumed together with the rotating wave approximation allows for the atomic state |1) to 
be coupled to one photon bath states and atomic state |0) to be coupled to bath states with 
two photons — the latter may be associated with the same mode or two different modes. 
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B. Coupled Amplitude Equations 



The time evolving state vector is written in the form: 

A 

+ EcoAAe-^*"^*|0)|2A)+ E coA;.e-^^'^^+'^'')*|0)|lAU, (4) 

A X,fj,,X<fi 

where C2,cia,coaa and cqa^j are interaction picture amphtudes for the various states. A 
straight-forward substitution into the Schrodinger equation leads to a set of coupled equa- 
tions for these amplitudes. 

The overall aim is to find the evolution of the atomic initial state, e.g. C2(t). We first 
obtain a set of coupled equations for the Laplace transforms of the amplitude equations, 



and then apply the approach of [33[ to make a convenient change of the amplitudes by 
incorporating the coupling constants via C2 = &2,cia = fi'A2&iA, cqa^j = 9x29,11^0x^1 (A < /x) 
and coAA = 9x29xiboxx- Laplace transforms are denoted bv a bar, and associated with the 
complex variable s. The coupled equations obtained are |4fil |: 

sb2{s)-l = -i'^gl2bix{s + i{ujx-uj2)) (5) 

A 

sbixis) = -i J2 9liboXfj,is + i{uj^-uJi)) 

/x,(^t>A) 

E 9liax,J}0fix{s + i{uj^ - uJi)) 

fj.,{lJ.<X) 

-igl^V2boxxis + i{ujx - ^^i)) - ihis + i{uj2 - ^^a)) (6) 

sboxx{s) = -iV2bixis + i{uji - ux)) (7) 

■s&oA/x(s) = -ibixis + i{uji - Uf,)) - iax,J>i,,{s + i{ui - ux)). (8) 

These involve a frequency independent parameter: 

9xi9^i2 /„x 

axfM = • (9) 

9X29,^1 

It is not possible to obtain explicit solutions for the new Laplace transform amplitudes. 
However, we can proceed further by first eliminating boXfi, boxx in the above equations, leading 
to a set of coupled equations for the remaining amplitudes b2{s) and bi^{s). 
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C. Integral Equation 



We then eliminate &2(s) in favour of (though not the reverse) giving a set of 

equations (see jl^) for the remaining amphtude bi^{s): 



i/s = is + i{uJx~uj2) + Y^ 

\ V 



9^1 



S + l[UJx + UJri — UJi — LJ2) 

2 



bix{s + i{ux - UJ2) 



V V 9X29,11 S + t{uJx + UJf,- iUi - iU2) s J 



(10) 



The above equation is for the amphtude 6ia(s) — it is of course coupled to all similar am- 
plitudes 61^ (s). 

By converting the sums to integrals, i.e. — ^ / duj^j_p{ijj^), we obtain a Fredholm 
integral equation of the second kind 



with 



f{uJx) + j diu^K{uJx,uj^) f{ujf,) = d{uJx) (11) 

/(cua) = hx{s + i{ujx-u:2)) (12) 

K{ujx,uj^) = B{uox,uo^)/A{uox) (13) 

d{ux) = CIA{ux) (14) 

A{ujx) = s + i{ujx- UJ2) + f dujr,p{ujr,) — T (15) 

J S + t[uJx + U^rf - UJi - UJ2) 

Biux,u,) = piu,) fe /""^^ . + ^) (16) 

\ 9X29fil S + t{uJx + UJ^-UJi-UJ2) s J 

C = -lis (17) 

The quantities f{u]x),K{u]x.,oj^^A{uJx),B{ijJx,ijJ,^,C and d{ujy) are all functions of the 
Laplace variable s. The quantity K performs the role of the kernel for the integral equation 
for the quantity /, and the quantity d makes the integral equation inhomogeneous. 

As pointed out in [4^, the integral equations for 6ia(s) and the other related equations 
for 62(5) and &oA^t(s) only involve the reservoir quantities gxi,9x2,P\ in terms of reservoir 
structure functions 

R,{ujx) = p{uJx)9li (18) 
i?2(^A) = p(^a)^7a2 (19) 



Since the atomic density operator 



PA = Trp 1^) (^1 



(20) 



also only involves reservoir structure functions and the reduced amplitudes b2{t),bix{t) and 
bo\fj,{t), an important result follows that the overall dynamics is entirely determined by the 
reservoir structure functions. 

The integral equation (fTTjl can be solved in different ways, such as: (a) numerical methods 
in which the integral equation is converted to a matrix equation (see next sections); (b) 
expansions using bi-orthogonal eigenfunctions (see Ref. |46|); (c) expansions such as the 
Fredholm expansion (see textbooks on integral equations). 

D. Numerical Solution of Integral Equation 

By discretising the mode frequencies uj\ the integral equation (lllj) can be converted to 
a matrix equation (the rows and columns are specified by the discrete frequencies and the 
integral over cj^ is approximated by a discrete sum) of the form 



which we then can solve for f . 

However, the obvious simple solution of the form f = (K + I)~^d is not correct. To get 
the time evolution of 62 for example, we need both its real and imaginary parts, b2r and 621- 
This means we need the separate Laplace transforms &2r(s) and &2i(s) to invert, and this in 
turn requires us to separately obtain fr{s) and fi{s). Thus we need the separate real and 
imaginary parts of /, both being functions of the complex Laplace variable s. These separate 
parts are obtained by analytic continuation of the real and imaginary parts of K{uJx,ujf^), 
f{uJx) and D{ujx) on the real s-axis. In other words, we first break the equation (K + I)f = d 
into its real and imaginary parts on the basis of the variable s being real, solve for fr{s) and 
fi{s) with real s, then obtain fr{s) and fi{s) for complex s by analytic continuation. 

In matrix form in terms of the real and imaginary parts, the integral equation ()21|) 
becomes 




(21) 




(22) 
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As emphasised above, the real, imaginary parts are denoted K^, Kj etc. and are identified 
for s real. Analytic continuation enables solution for and fj which will apply for all s. 

Solving for and fj by matrix inversion, gives the LT amplitude bi\{s + i{ujx — uj2))- 
The coupled amplitude equations © give 62(5). For the discretised frequency form, the real 
and imaginary parts can be written as the scalar products 

b2r{s) = (l + r-^)/s (23) 

b2^is) = -V-Tr/S, (24) 

where r = {pAfi'A2}- The amplitude b2{t) then follows from a numerical inverse Laplace 
transform of both b2r{s) and b2i{s). 

III. APPLICATIONS 

A. Lorentzian resonances 

Before looking at two special cases, we will first generate some results for a reservoir 
structure composed of a sum of Lorentzian resonances. The resonances need not be well 
separated from each other in frequency space, but must be well above zero frequency. The 
general form of structure function is then: 

pgliu) = ,pgliu) = E ^ ■ (,_,c).V(r„/2)^' (25) 
where a simple scaling 77 has been introduced between the strength of coupling for the upper 
and lower transitions, gl^ = ?75'a2- this situation ax^j, = 1, [see equation Q]. The ath 
Lorentzian is weighted by Za (not necessarily positive) and has a resonance centred at the 
frequency with width Fq,. 

For these cases of reservoir structure function, the integral in equation ()15|) is straight- 
forward and gives for the function A 

AM = s + i{ux - ^2) + E \ c ^ ■ (26) 

In what follows we define the reservoir frequency and resonances relative to the average 
of the transition frequencies, i.e. we let 

AUO = UOx- {OJI + Co'2)/2 

Acu' = uj^-{loi + uj2)/2 (27) 



and define 

5^ = ujI- (c^i + uj2)/2. (28) 

The detuning 5a is then the offset of the ath Lorentzian from the average frequency such 
that if 5a = ±5 the ath Lorentzian is resonant with tui, or uj2i according to the sign of 5. 
Then, together with equation for 5, we find in this notation that 



A(Aa;) = s + z(Au; + 5) + 



^ s + r„/2 + i(Aw + 5„) 



(29) 



and 



1 1 

r]s s + z(Aco' + Au;') 



27r (Ac^'-5,)2 + (r,/2)2" 



(30) 



As before, the kernel is then /^(Acj, Au;') = S(Acj, Au;')M(Acj). 



B. Case A: High-Q Cavity Model 

In this section we consider the special case of the three-level atom coupled to a single 
high-Q cavity mode. For simplicity, identical coupling constants will be assumed. The atomic 
transition frequencies may be unequal. Thus we have gxi = g\2 = g\, and r/ = 1 in equation 
(|25|l . In this case the frequency dependence of the reservoir structure function will be due 
to resonant behaviour of the coupling constants, as a quasimode theory of such systems 
demonstrates Thus, parameterising the single weight Zi by Q"^ we have: 

Ri = R2 = R = pxgl (31) 

_ m2 1 

~ ^ ■ (Aoj-(5)2 + (r/2)2' ^^^^ 

where we denote the single detuning 5i by 5. 
The kernel K [see equation (fT!?|l ] is given by: 

K(A A ') + ^^^^ + '^^ + ^/^^^^^ + ^(^^ + ^^'^ 

^ ' ^ 27r s((Acu'-5)2 + (r/2)2)(s + ^(Acu + Acu'))g(Acu)' ^ ^ 

where (5(Acij) = [s + z(Au; + 5)][s + z(Aci; + 5) +r/2] +fi2. The frequency integral in equation 
fll3|) is integrated over the range ±oo for convenience. The matrices K^, Kj are determined 
from the kernel K and the matrices d^, dj from an analogous equation for d. The amplitude 
62 (t) is obtained by numerical integration. 
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FIG. 2: High-Q cavity model. Time evolution of the probability of finding the system in state 2; 
P{t) = |62(i)P- The reservoir structure function is given by equation H32|) . with parameters F =1 
and: (a) ^1 = 5.0; (b) 17 = 1.0; and (c) 0, = 0.5, in scaled units. The grid size for the discretised 
kernel was 150x150 chosen with a range of ±30 for Auj and Auj' in scaled units. The case of 
resonance is shown, where each transition is resonant with the reservoir structure: S = 6 = 0. 
Results were computed with the integral equation method of section Hi Dl 

We note that in the special case of a structure resonant with both the transitions, ui = 
002 = ^0 = and 6 = 6 = 0. The kernel then reduces to the case studied in |46|, i.e. 




where now Q{Auj) = {s + iAuj){s + iAuj + T/2) + fi^. 

Typical results for the time-dependent upper state probability are shown in figure 12 
where the parameters used are F =1 and (a) Q =5 (b) Q =1 (c) Q =0.5 (effectively in scaled 
units). The damped oscillation in case (a) is typical of non-Markovian decay. In case (b) 
the oscillations weaken. Further reduction in coupling in case (c) removes the oscillations 
and the decay is closer to exponential. 

The case just described applies when both transitions are coupled to a single reservoir. 
A different situation applies when the two transitions are coupled to two separate reservoirs, 
and the results for this case are presented in The equations are then simpler, and the 
integral equation for the amplitudes can be solved analytically. The decay of the excited 
state is now less oscillatory. Unlike the single reservoir case, the two photons emitted can 
now be distinguished. 

The typical effects of detunings on the upper state probability are presented in figures El 





Ffi^ (s + iAuj + +F/2) (2s + i{Au + Acj') 
2ti s(Aw'2 + (F/2)2)(s + z(Aa; + Auj'))Q{Auj) ' 



(34) 
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FIG. 3: High-Q cavity model. Time evolution of the probability of finding the system in state 2; 
P{t) = |62(t)p. The reservoir structure function is given by equation H32|) . with parameters F =1 
and 0, = 5.0 in scaled units. The atomic transitions have equal frequency: 5 = 0. The centre of the 
resonance is offset from the atomic frequency by the detunings in (a)-(d) such that 5 = 0,5, 10, 20. 
The results for (a) and (b) are computed using the integral equation method of section IIIDI and 
the results for (c) and (d) are computed using the master equation method discussed in section HVl 

and El Figure IHl shows the case where the two atomic transitions have identical frequency 
and are detuned from the reservoir structure by increasing amounts [curves (a)-(d)]. It is 
seen that the oscillation frequency initially decreases, but then increases with the detuning. 
The damping of the oscillations remains about the same, but the decay of the population 
itself reduces as the detuning increases. For the large detuning cases 6 = 10, 20 (curves (c), 
(d) in figure El) the integral equation method was very slow, and therefore these curves were 
computed using the master equation method to be discussed in Section HVl As we will see, 
this method is equivalent to the integral equation method since the same reservoir structure 
function is involved. Figure shows the situation for a fixed cavity mode detuning 6 = 4, 
but with atomic transition frequencies that may be different. The chained and dotted curves 
in figure m show the effect of unequal atomic transition frequencies. In fact, since 6 = ±5 
for these two cases, we can see from equation (jT}, and 6 used in equation (j32ll . that 6 = ±6 
corresponds to the Lorentzian structure being resonant with the lower atomic transition, 
or the upper transition, respectively. The effect on the population dynamics is similar in 
both cases — the population is lost slightly less quickly, because of less efficient coupling, and 
the frequency of the oscillations changes. The frequency of population oscillation is reduced 
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FIG. 4: High-Q cavity model. Time evolution of the probability of finding the system in state 2; 
P{t) = |62(i)P- The reservoir structure function is given by equation H32|) . with parameters F =1 
and = 5.0 in scaled units. The detuning (5 = 4 and the off-set is (5 = (solid line), +6 (chain), 
and —5 (dotted). The results are computed using the integral equation method of section Hi D[ 

when the structure is resonant with the lower transition, and increased when it is resonant 
with the upper transition. Note that for equal atomic transition frequencies both figure El 
and figure m display a case where the population drops close to zero periodically. This is for 
5 = 5,4 respectively. 



C. Case B: Photonic Band-Gap Model 

In this section we consider a simple model of a three-level atom coupled to a photonic 
band-gap system. For simplicity, we again assume equal coupling constants and transition 
frequencies. Thus gw = g\2 = gx, so that t] = 1. In this case the rapid frequency dependence 
of the reservoir structure function is due to mode frequency gaps, as a quasimode theory of 
such systems demonstrates single reservoir structure function given as the difference 

of two Lorentzians is assumed ;34j to approximately represent the frequency gap, which is 
assumed to be detuned from the average atomic transition frequency by 6 (which is the same 
for both Lorentzians). In this case we let Zi = and Z2 = —^2 with 00^ = 002 = 00'^, so 
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that 6i = 62 = S. Then the reservoir structure function is given by: 



Ri = R2 = R = pxgl (35) 



2ti {ux - uj-)^ + (ri/2)2 27r (cux - uj-f + (12/2)2 ' 

P| r2 



(36) 



27r (Acj - 5)2 + (ri/2)2 271 (Acj - + (13/2)2 ' 
The Lorentzians have couphng strengths i^i, ^2 and reservoir structure widths Fi, r2, where 

f = f (37) 
J- 1 J- 2 

< Ti. (38) 

These conditions are required so that R is always pos itive and has a zero where Acj = 5. This 
locahsed zero models the photonic band-gap fsee IsJ-IsS^). The atomic transition frequencies 
might not be equal to each other, with the difference being represented by the parameter 5, 
equation (0), but if 5 = 5 = then both transition frequencies would be resonant with the 
zero in the reservoir structure function equation (j36|) . 

The kernel K and the function d are now given by (see equations ()13fl4j) ): 

1 iSllTi - nlVi) (Acu' - 5f 



271 ((Aa;' - Sy + (ri/2)2)((Aa;' - 5)2 + (12/2)2) 
{2s + i{Auj + Au')) 
s{s + i{Auj + Au')) 



X 



{s + zjAcu + S) + r,/2)is + zjAuj + 6) + T2/2) 

Qs{s + iiAu + 5)) ^ ^ 

. z!(- + -(^- + ^) + r./2)(- + .(A. + 5) + r2/2) 

^ ^ s Qsis + t{Auj + 6)) ^ ^ 



with 



Qsix) = [x + ^(5-(5)](x + ^l/2)(x + ^2/2) + x(^]?-^]^). (41) 

In equations (jHIIjl - PTj) we have used relation (jSTj). 

The matrices K^, Kj are determined from the kernel K and the matrices d^, dj from an 
analogous equation for d. The amplitude b2(t) is obtained by numerical integration. 

Typical results are shown in figureElfor the case of two equal atomic transition frequencies, 
also coincident with the gap frequency. Thus ui = uj2 = uj'^, so that 6 = 6 = 0. The 
parameters used are Fi = 4 and Qi = 1 with r2 =1.0, 1.5 and 2.0. Non-Markovian decay is 
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FIG. 5: Photonic band-gap model. Time evolution of the probability of finding the system in state 
2; P{t) = |62(i)P- The reservoir structure function is given by equation ()36() . with parameters Fi 
=4 and = 1 and: (a) F2 =1 (b) F2 =1.5 (c) F2 =2.0 in scaled units. The case of resonance is 
shown, where each transition is resonant with the gap: 6 = 5 = 0. 

seen for small t as the initial quadratic behaviour. Population trapping effects are found in 
all three cases. The amount of trapped population increases as r2 increases. 

Figure El shows the effect of detuning the two (equal) atom transition frequencies from 
the centre of the gap. The parameters used are as in figure El with F2 =1. Again, for the 
large detuning case 6 = 8 (curve (e) in figure E)) the results were calculated using the master 
equation method of Section HVl We see that there is a loss of the population trapping effect, 
although for larger detunings 6 the rate of decay is considerably reduced. However, this 
reduction can be attributed to the decrease in the density of states when detuned far from 
the gap and outside the scope of the second, positive, Lorentzian. In that sense it is an 
artifact of the model. The effect is eliminated for the parameters shown in figure [7| where 
Fi is comparatively large. The parameters are as in figure El except that Fi = 50 and the 
detunings selected are different. All the graphs in figure El were computed using the master 
equation method. In this case only the loss of the trapping is seen as the detuning increases, 
and the curves appear to saturate as the background density of states flattens out for the 
larger detunings. 
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FIG. 6: Photonic band-gap model. Time evolution of the probability of finding the system in state 
2; P2{t) = |&2(i)P- The reservoir structm'e function is given by equation H36|) . with parameters Fi 
=4 and J7i = 1 and F2 =1 in scaled units. The atomic transitions have equal frequency: 5 = 0. The 
centre of the gap is offset from the atomic frequency by the detunings in (a)-(e) of 5 = 0,0.5, 1,4 
and 8. The results for (5 = to 4 are computed using the integral equation method of section Hi Dl 
and the result for (5 = 8 is computed using the master equation method discussed in section ITVl 




10 20 30 40 50 

t 



FIG. 7: Photonic band-gap model. Time evolution of the probability of finding the system in state 
2; P2{'t) = \b2{t)\'^. The reservoir structure function is given by equation (|36() . with parameters 
Fi =50 and Jli = 1 and F2 =1 in scaled units. The atomic transitions have equal frequency: 
6 = 0. The centre of the gap is offset from the atomic frequency by the detunings in (a)-(e) of 
5 = 0, 0.25, 0.5, 1 and 2. These results are computed using the master equation method discussed 
in section HVl 
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IV. MASTER EQUATION THEORY 



A. Pseudomode Treatment 



We have previously found that Lorentzian reservoir structures can be connected, in a non 



Mia 



perturbative way, to Lindblad master equations using the pseudomode approach 
which involves considering the positions and residues of the poles of the reservoir structure 
function in the complex frequency plane. Each pole is associated with a pseudomode. 

This approach was applied to the case of a two-level atom in a high-Q cavity with a 
single resonance of width Fi The system takes on the dynamics of a damped 

Jaynes-Cummings model where the atom exchanges energy with a pseudomode, from which 
the energy is lost. In the case of our three- level atom in a high-Q cavity the pseudomode 
approach would result in a master equation of the form 



dtp 



ttCAV - 
^0 



(^a[aif) — 2aif)a\ + f)a\aij , 



(42) 



with the Hamiltonian 



CAV 




+ (cji + cj2)|2)(2| + uj'a\ai + H'x'f 



(43) 



Here the term H^^l' represents the atom-pseudomode coupling 



ffi'^ = /ifi(al|0)(l| + ai|l)(0| +al|l)(2| + ai|2)(l 



(44) 



Comparison of the calculated results shows that this master equation reproduces the results 
of Section Fill Bl for the atomic populations obtained by the integral equation method. 

The pseudomode approach was also applied to the case of a two-level atom in a photonic 
band-gap system [24 . In the case of cascade atom in the band-gap material, as in Section 
nil (^1 this approach 2^, 0] would, in general, result in a two-mode master equation of the 
form: 



dt 



P 



n 

T' 
~2 



^0 >P 



£1 

2 



2aif)a\ + f)a[ai 



— f- (4«2P - 2a2pal + pala-^ 



(45) 



17 



with the Hamiltonian 



+ (cui + co'2)|2)(2| + u^alai + ^"02(12 



+fi — - — {a[a2 + aifla 
+hSI- (46) 

This time the atom-pseudomode couphng is to the pseudomode labelled '2', i.e. 

= ^^PBG (4|0)(1| + a2|l)(0| + a^|l)(2| + a2|2)(l|) , (47) 

similarly to equation (jlH). The two modes (or pseudo modes) of this case are coupled to 
each other in the Hamiltonian For the case when the reservoir structure function has a 
zero, as has been prescribed by equation (j?fj) . we find that F'^ = 0. However, Fg is non-zero 
and is related to the large u behaviour of the reservoir structure function R. This parameter, 
and the atom-pseudomode coupling strength are given by 

r'2 = Ti + F2 , (48) 



^PBG = V^i - ^2 (49) 

which are both positive. Again, comparison of the calculated results shows that the master 
equation reproduces the results of Section UlI (^I for the atomic populations obtained by 
the integral equation method. For the full density matrix, 36 coupled differential equations 
are required with the two excitations produced from the three-level ladder system. 



B. Quasimode Treatment 

Markovian master equations describing Lorentzian reservoir structures can also be ob- 



tained based on quasimodes 



4^ . which are identified with approximate modes of the 



EM field. The basic concepts and equations for quasimode theory are set out in the Ap- 
pendix. For the present paper we consider the case where the continuum quasimode density 
Qci^) and the n discrete-continuum coupling constants Wi{A) are slowly varying functions 
of the continuum quasimode frequency A. In this situation we obtain Markovian master 
equations (equation (jH21)) for the density operator p describing the system consisting of the 
atom plus the n discrete quasimodes. The continuum quasimodes constitute the reservoir. 
In addition simple expressions are obtained for the reservoir structure function Rk{oj) for 
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the true modes [equation (|HH)]. These results will now be applied to the high-Q cavity and 
photonic band-gap cases. 

For the high-Q cavity case, only a single discrete quasimode is involved. The master 
equation for the cascade atom plus single discrete quasimode system is then 

d —i-^ r 

= y [Hs,p\ + 2 {[^I'M] + [aip,d\]}, (50) 
where the decay rate V is given by 

V = 2t,q,\Wi\\ (51) 

and the system Hamiltonian is 

Hs = J2 Vk fiuJk {at^k - ^k^k) + a! 

k=l,2 

+ y: in\i,d,at + H.c), (52) 

k=l,2 

with at = |1)(0|, = |2)(1|, ^71 = § + V2 = ^ + The atomic Hamiltonian is that 
given in equation (j21), apart from an additive constant {2huji + huj2)/3. The choice An = 
A21 = ^ (real) is natural since we are dealing with a situation where both transitions have 
equal coupling constants with the EM field, and we choose z/i = u'^ to have the quasimode 
frequency coincide with the peak in the reservoir structure function. These choices then 
result in the same master equation as equation (j42|l . A straightforward evaluation of the 
reservoir structure function from equation ()84j) gives 



for both transitions. This result is the same as in equations ()3ipH2|) . Thus the high-Q cavity 
case can be treated via quasimode theory involving a single discrete quasimode, which may 
be identified as the cavity mode. 

For the photonic band-gap case, two discrete quasimodes are involved. There is now 
more choice for the parameters in the quasimode model and it turns out we can choose 



these to obtain the same reservoir structure function as in equations ()35p36|) as well as 
generating a master equation which is the same as equation obtained from pseudomode 
theory. However, to achieve this we must first introduce two new discrete quasimodes, whose 
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annihilation operators bi, 62 are defined by 

b^ = n~y^Trp^(-W2ai + Wid2) (54) 
$2 = K'^^'^TTp^{W*ai + W;d2) (55) 

Av = i(ri + r2), (56) 

where the decay rates are 

Ti = 27Tg, \Wif (57) 

T2 = 27rg, \W2f. (58) 

Note that the 61, 62 and their adjoints satisfy the standard Bose commutation rules, which 
is a consequence of the matrix involved in (jHH l55|l being unitary. 

From equation ()82j) we immediately see that the master equation is given by 

j/ = Y + 1^^' + mKpbi] + [hpM, (59) 

showing that relaxation only involves the new quasimode with annihilation operator 62? the 
decay rate being the sum of the original rates Fi and r2. This result for the relaxation terms 
is the same as in the pseudomode master equation P5|) . 

The system Hamiltonian Hs can then be written in terms of the new quasimode operators 
by inverting equations (jnHES)) to give di, 0,2 in terms of 61, 62 and substituting in equation 
fl83|) . There are still many choices that can be made for the quasimode frequencies Ui 
{i = 1,2), the magnitude and phases of the atom-quasimode coupling constants Xki {k = 
l,2;i = 1,2), the quasimode-quasimode coupling constant Vu and the phase for discrete- 
continuum coupling constants Wi {i = 1,2). The magnitude of the Wi are equivalent to the 
decay rates Fj. The atomic term in the system Hamiltonian Hs is unaltered 

^51 = J2 Vk fi^k - ?fe )• (60) 

fc=l,2 

Since the photonic band-gap model involves two Lorentzians with the same centre fre- 
quency u'^ it is natural to choose the discrete quasimode frequencies Ui to be both equal 
to uj'^, and since both atomic transitions have equal coupling constants with the EM field 
it is natural to choose the atom-quasimode coupling constants to be independent of the 
transition. Thus 

I'l = 1^2 = (61) 
A/ci = Ai, Ajt2 = A2. (62) 
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We also introduce the various phase factors via the definitions 



Vi2 = \Vi2\ exp(i^i2) 

Wi = \Wi\ exp(i0i), W2 = \W2\ exp(z02) 
Ai = |Ai| exp(i6'i), A2 = IA2I exp(i02)- 



(63) 
(64) 
(65) 



In view of fl61|) and the unitary relation between the new and original discrete quasimode 
operators, it follows that the quasimode energy term in Hs is given by 



Hs2 = hu'{b{h + bl 62). 



To obtain the required atom-quasimode interaction we choose 



\Vl2\ 
^12 



2vrir2 



TT 



ii-02 + (2z/ + l)- (z/ = 0,±l,±2,..) 



TT 



(/.1 + 02 = (2z/+l)--/i27r (/i = O,±l,±2,.0 
and find that the quasimode- quasimode interaction is given by 

Hs3 = fil^/rj\{b\b2 + blh). 



(66) 

(67) 
(68) 
(69) 

(70) 



To obtain the required atom-quasimode coupling together with condition ()62j) we choose 

IA2I =^^2. 



lAii = nu 



Ol- (pl = 02- 02 



2 ' 
2 



1 



i2Tx (e = 0,±l,±2,..) 



(71) 
(72) 
(73) 



where the (real) quantities Vti, VL2, Fiand r2 are related as in equations (jSHISHl)- We thus 
find that 

Hsi = n QpBG {h {^t + (^t)+ bl (?r + ^2")) (74) 
where, as in equation (jlUjl . 

' (75) 



^PBG = V ^1 — ^2- 

Overall, with Hs = Hsi + Hs2 + Hs3 + Hsa the system Hamiltonian (and hence the master 
equation) is the same as that obtained via the pseudomode approach. The key feature is 
that we started with two original coupled discrete quasimodes with equal frequencies, whose 
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parameters satisfy the conditions in equations IHHl IHZl EH EHl IZH IZ21 ESI)- Two new- 
discrete quasimodes were introduced, also coupled and with equal frequencies. However, the 
first of these new quasimodes bi is neither coupled to the atomic transitions nor involved in 
relaxation. The other new quasimode 62 is coupled to the atomic transitions and is involved 
in relaxation with a rate equal to the sum of the original rates. 

Finally, with the conditions referred to in the last paragraph applying, it is a straightfor- 
ward process to evaluate the reservoir structure function from equation ()84|). We find that 
this is the same for both transitions and is given by 



Rkioj) = R{uj 



^Ti 1 



27r (cu - cj'=)2 + (ri/2)2 

""^ 2^ {uj-u^Y + {V,/2f ^^^^ 
Thus the photonic band-gap case can be treated via quasimode theory involving two coupled 
discrete quasimodes of equal frequency, where one may be identified as a band-gap mode 
and the other background mode. 



V. SUMMARY 



The dynamical behaviour of a three level cascade atom coupled to a structured reservoir 
(typically of EM field modes), and initially in the upper state has been analysed via Laplace 
transform and exact master equation methods. In the Laplace transform approach, the 
atomic density operator is determined from the solutions of integral equations, in which the 
properties of the structured reservoir only appears via reservoir structure functions, all es- 
sentially given by the product of the mode density times the square of coupling constants. In 
the cascade system two distinct reservoir structure functions are involved. The dependence 
of the dynamics solely on reservoir structure functions is required for treating structured 
reservoir problems via pseudomode theory, so our results suggest that it may be possible to 
extend pseudomode theory to problems involving more than a single photon excitation of 
the reservoir. 

Solutions of the integral equation involve discretising the frequency space into points, 
so the matrix inversion step involves a matrix with AN"^ elements compared to, say, 0{N^) 
elements represented by the original coupled amplitude equations in an equivalent discretised 
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form. Considerations based on analytic continuation were required for obtaining solutions 
of the problem using numerical methods, including carrying out an inverse Laplace trans- 
form. For some oscillatory problems (specifically, with large detunings) there were numerical 
difficulties with the integral equation approach. 

Two important physical situations have been studied using these numerical methods. 
The first applies to a simple model of a three-level atom with identical cascade transitions 
coupled to a single high-Q cavity mode, involving a single Lorentzian reservoir structure 
function. Non-Markovian oscillatory decay of the excited state has been demonstrated. The 
second applies to a model of a three-level atom with identical cascade transitions coupled to 
a photonic band-gap system, involving a single reservoir structure function modelled simply 
as the difference between two Lorentzian functions. Again we have shown the non-Markovian 
oscillatory decay of the excited state and population trapping effects. 

We have also studied both the high-Q cavity and the photonic band-gap cases via the 
master equation method. The master equation can be obtained both from the pseudomode 
approach and the quasimode approach, and numerical calculations of the upper state proba- 
bility via the master equation method give the same results as the integral equation method 
based on essential states. This shows that complicated non-Markovian decays into struc- 
tured EM field reservoirs can be described by Markovian models in which the atomic system 
is augmented by one or two pseudomodes or quasimodes, which in the quasimode approach 
themselves undergo Markovian relaxation into a fiat reservoir. For the high-Q cavity case 
the single pseudomode or quasimode may be identified as the cavity mode, for the photonic 
band-gap case the two coupled pseudomodes or quasimodes may be identified as band-gap 
and background modes. 
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APPENDIX : QUASIMODE THEORY 



In quasimode theory |43|, an atomic system with upward, downward transition opera- 
tors , and transition frequencies Uk is coupled to a set of n discrete bosonic quasimodes 
with annihilation, creation operators a,, a| and frequencies Ui. The coupling constants are 
denoted Xki. The discrete quasimodes may be coupled to each other with coupling constants 
Vij. In addition, there is a coupling between the discrete quasimodes and a continuum of 
quasimodes whose annihilation, creation operators are 6(A), b{Ay and whose frequencies 
are A. All coupling terms are in the RWA. The continuum quasimodes have mode density 
f)c(A) and the coupling constant with the discrete quasimodes are denoted Wi{A). The 
continuum quasimodes are neither coupled to each other nor to the atomic system. In Refs. 



44l | the quasimodes were taken to be approximations to the true EM field modes. The 
overall Hamiltonian is given by 

H = Ha + Hq + Haq, (77) 

with 

Ha = fiujk {Bt^k - K^k) (78) 

k 

i i+i 

+ T.Jd^ Qc{A)[nW,{A) a\ 6(A) + H.c] 

i 

+ J dAg,{A)nAb{Ayb{A) (79) 

HaQ = Y.Y.i^K^^^^t + H.c). (80) 

k i 

In equation (f75j) the quantities 7]^ are chosen so that Ha equals the atomic Hamiltonian 
given in equation apart from an additive constant {2huJi + huj2)/3. The annihilation 
and creation operators for the discrete quasimodes satisfy the Kronecker delta commutation 
rules, whilst those for the continuum quasimodes satisfy Dirac delta function commutation 
rules: 

[ai,a]] = 6ij 

[6(A),6(A')t] = <5(A-A')MA). (81) 

The Qc factor on the right side of the last equation gives annihilation and creation operators 
that are dimensionless. 
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In Refs. the true EM field modes were described by annihilation and creation op- 

erators A{uj) and A{ujy that were obtained by Fano diagonalisation methods. This involved 
writing A{uj) as a linear combination of the discrete quasimode annihilation operators di plus 
an integral of the continuum annihilation operators b{A), and requiring that the true mode 
annihilation operator is an eigenoperator for the quasimode Hamiltonian ()79|1 with energy 
hu. Expressions for the atom-quasimode coupling ()80|) could then be written in terms of 
the true mode annihilation and creation operators (equations (18), (19) of Ref. 4J|) and 
a result for the reservoir structure function Rk{uj) associated with the k atomic transition 
obtained (equation (50) of Ref. j^|). These results applied to the general case where the 
continuum quasimode density gd^) and the discrete-continuum coupling constants Wi{A) 
were not necessarily slowly varying functions of A. 

For the present paper we are interested in the situation where the dynamics of the system 
consisting of the atom plus the n discrete quasimodes is Markovian. The reservoir consists 
of the continuum quasimodes. Markovian behaviour occurs when gdA) and the VFi(A) are 
slowly varying functions of A and can be treated as constant. The flatness of the continuum 
quasimode density and the coupling constants results in the reservoir correlation time being 
short enough for Markovian behaviour to occur. The density operator p for the system 
consisting of the atom plus the n discrete quasimodes can then be shown to satisfy the 
Markovian master equation 
d ^ —i 

where the system Hamiltonian is 



-p = — [Hs, p]+Yl Wi W; { [a„ paj] + [a,p, aj] }, (82) 



Hs = Vk fiuJk [Bt^k - ^k ^k) + ^1 



Y m, d\ a, + E E (^K^ai Bt + H.c). (83) 

i^i k i 



In this case of a fiat continuum the reservoir structure function is given by 

Rk{uj) = gc 



|2 ' 



(84) 



where 



(85) 
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Pnico) = \ujB^ -n\- iTiQ, W, (u;E„ - W;, (86) 

and the n x n matrix ft has elements 

nij = UiSij + {l-Sij)V,i, (87) 
with the adjugate matrix (cuE^ — O)"^-^"^ related to the inverse via 

n obtaining this result we note that the frequency shift matrix Fij{uj) (equation (27) of Ref. 
4^ ) is zero. 
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